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ABSTRACT 

We present an implementation of the Schwarzschild orbit superposition method which can be 
used for constructing self-consistent equilibrium models of barred or non-barred disc galax¬ 
ies, or of elliptical galaxies with figure rotation. This is a further development of the publicly 
available code SMILE; its main improvements include a new efficient representation of an 
arbitrary gravitational potential using two-dimensional spline interpolation of Fourier coef¬ 
ficients in the meridional plane, as well as the ability to deal with rotation of the density 
profile and with multicomponent mass models. We compare several published methods for 
constructing composite axisymmetric disc-bulge-halo models and demonstrate that our code 
produces the models that are closest to equilibrium. We also apply it to create models of 
triaxial elliptical galaxies with cuspy density profiles and figure rotation, and find that such 
models can be found and are stable over many dynamical times in a wide range of pattern 
speeds and angular momenta, covering both slow- and fast-rotator classes. We then attempt to 
create models of strongly barred disc galaxies, using an analytic three-component potential, 
and find that it is not possible to make a stable dynamically self-consistent model for this 
density profile. Finally, we take snapshots of two A^-body simulations of barred disc galaxies 
embedded in nearly-spherical haloes, and construct equilibrium models using only informa¬ 
tion on the density profile of the snapshots. We demonstrate that such reconstructed models 
are in near-stationary state, in contrast with the original iV-body simulations, one of which 
displayed significant secular evolution. 
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1 INTRODUCTION 

An important task that often arises in stellar dynamics is the con¬ 
struction of self-consistent equilibrium models of galaxies that have 
a given density profile p{x) and, optionally, satisfy certain kine¬ 
matic constraints. Self-consistency means that the gravitational po¬ 
tential is determined by the distribution of the mass in the galaxy 
via the Poisson equation. Equilibrium models satisfy the collision¬ 
less Boltzmann equation with a steady-state distribution function. 
In practice, we usually require the models to evolve very slowly 
compared to the dynamical (crossing) timescale. 

There exist numerous methods for creating (near-)equihbrium 
galaxy models of varying complexity and geometry. After re¬ 
viewing of these approaches in Section with a particular fo¬ 
cus on the methods suitable for disc galaxies, we concentrate on 
the Schwarzschild orbit superposition method. Up to now, it was 
mostly applied to elliptical galaxies or to the inner parts of the 
Milky Way. We describe the challenges that need to be met in its ex- 
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tension to highly flattened, non-axisymmetric stellar systems with 
figure rotation, and present our implementation in Section]^ Then 
in Section |4] we apply it to construct models of some analytical 
density profiles that are frequently used in the literature, in particu¬ 
lar, comparing several published methods for creating equilibrium 
models of axisymmetric disc galaxies, and create models of barred 
disc galaxies based on the density distributions taken from W-body 
simulations. We discuss our results and possible future applications 
in Section|5] and summarize in Section|^ 


2 METHODS FOR CREATING EQUILIBRIUM MODELS 

The focus of this paper is on disc galaxies, therefore we restrict our 
review to the methods suitable for such systems which are typically 
multicomponent (include a strongly flattened disc, which may or 
may not be axisymmetric, a rather spherical halo component, and a 
bulge or a bar that might well be triaxial). 

Most of the methods specifically tailored for disc galaxies pro¬ 
duce axisymmetric models that are close but not exactly in equi- 


© 2015 RAS 


Schwarzschild method for disc galaxies 2843 


librium. iHemguisl ( Il993h introduced a method based on the lo¬ 
cally Maxwellian approximation to the velocity distrihution func¬ 
tion and a prescribed density distribution for each component. It 
used the spherical (for the halo) and axisymmetric (for the bulge) 
Jeans equations and an approximation of constant radial-to-vertical 
velocity dispersion ratio for the disc together with epicyclic approx¬ 
imation for the azimuthal component of velocity dispersion. In this 
method, all component s are spherical except, in some part s of the 
analysis, the disc; later lBoilv. Kroupa & Penarrubial l l200lh gener¬ 
alized this method to allow for axisymmetric halo and bulge, while 
using an approximate method to squeeze the shape of velocity el¬ 
lipsoid proportionally to the flattening of the potential. A drawback 
of this approach is the use of Maxwellian velocity distrihution func¬ 
tion, which may not always he appropriate fo r an equilibrium solu¬ 
tion teazantzidis, Magorrian & Moorell^004h . 

iKuiiken & Dubin^ jl995ti used a different approach, starting 
from a suitable anzatz for the distrihution functions of each compo¬ 
nent, which used two classical integrals - energy and 2 -component 
of angular momentum for the spheroidal components, and added 
the third approximate integral (energy of vertical motion) for the 
disc. The density profile corresponding to the given distribution 
function was then computed iteratively, taking into accou nt the po¬ 
tentials of all components. IWidrow & Dubin^ ( l2005h modified 
this scheme to allow for a cuspy density profile and a central mas¬ 
sive black hole, while retaining analytic expressions for the distri¬ 
hution functions. However, with this approach the density profile 
of the halo and bulge are determined im plicitly and are hard to pre¬ 
scribe ; an improvement introduced by IWidrow. Pvm & Dubins^ 
1 I 2 OO 8 I) instead derived the isotropic distribution function for the 
given density profile numerically, using the Eddington inversion 
formula in a spherically-symmetric approximation of the total po¬ 
tential. A similar approach (iterative solution for the potential 
for given analytically described distribution functions for each 
component) is employed in the widely used Besangon model of 
Milky Way ( IRobin et al.l |200^ . Binney and collaborators advo- 
cate the use of actions as arguments of distribution function (e.g. 
iBinnev & McMillanl201 ll) ; however, their models are not designed 
to be self-consistent dynamically and instead focus on the fitting of 
kinematic properties of observed stellar populations in the Milky 
Way disc. 

Another possible approach is to build up the composite sys¬ 
tem by introducing the components of gravitational potential one 
by one, i n the adiabatic approximation. In the first such implemen¬ 
tation byi arne i{I^, one begins with W-body realizations of 
two separate analytic spherical equilibrium models for bulge and 
halo, superimposing them on top of each other and allowing them 
to relax. This first stage, however, induces non-negligible changes 
to the original density profiles as they reach the new equilibrium in 
the joint gravitational potential. Then the analytically defined disc 
potential is grown adiabatically, without actually adding any disc 
particles, but allowing the particles in the spheroidal components 
to adjust to the combined potential. Finally, the disc is populated 
with particles using the loc ally Maxwellian distribution function. 
iMcMillan & PehnM ( |2007|) improve this method by constructing 
the spheroidal components from an exact distribution function de¬ 
rived from an anisotropic generalization of Eddington inversion 
formula, using the monopole terms of potential of all three com¬ 
ponents. Then only the non-spherical part of the disc potential is 
adiabatically introduced into the live W-body system of bulge plus 
halo, and finally the disc particles are populated from a more ac¬ 
curate distribution function, avoiding the Maxwellian approxima¬ 
tion. This method offers little control over the shape of the halo 


or bulge (which slightly flatten during adiabatic contraction), but 
yields close-to-equilibrium configurations even for rather “warm” 
discs. 

All of the above mentioned methods are designed for ax¬ 
isymmetric potentials only, often make an implicit assumption 
that the motion of stars in the disc com ponent is integrable 
(e.g. in the torus construction method of iBinnev & McMillanI 
l201lh . and usually have a rather restricted choice for the mass 
and velocity profiles of various components. It has long been 
recognized that non-axisymmetric features in galactic discs, 
such as bars, spiral arms and rings, give rise to a variety 
of p henomena asso ci ated with re sonances and chaotic orbit s 
(e.g. IPehnenl 120001 ; IFuxI 1200 ll ; iRomero-Gomez et alj l2006l) . 
Many studies have been devoted to the detailed analysis of 
orbital structure of ba rred disc galaxies described by analyt¬ 
ical potentials (e.g. I Contopoulos & Pap avann opoul osl [l980l; 


Athanassoula et al. 198^ 


.,onto£oulos & Pap avann oDoul osI I19o0l; 
I; LPfenn igein i984j; iTeuben & Sand^ 


19851 ; ISkokos. Patsis & Athanasso ula l2002l) or W-body mod- 


els (e.g. IPfenni ger & Fri e^ 119911; iBeren tzen et a lT jl998l ; 


iHarsoula & Kalapotharakos 20091 ; Fragkoudi et all 20151) . To 

address the question which of the many kinds of orbits actually 
contribute to the self-consistent potential, so-called response 
models are constructed , follo wing the method developed by 
IContopoulos & Grosbdll (Il988l) . It starts from an assumption for 
the gravitational potential (including non-axisymmetric perturba¬ 
tions and a value for pattern speed), then the periodic orbits in 
the given potential are computed and superimposed to obtain the 
“response density” profile. The parameters of the potential are then 
iteratively adjusted until the response density matches the poten¬ 
tial. A number of indicators are used to determine self-consistency, 
including the amplitude and phase of various angular harmonics as 
functions of radi us. This technique has been appli ed to barred spi¬ 
ral galaxies (e. g. iKaufmann & Contopoulosll 19961); a mo dification 
introduced in Kalapotharakos. Patsis & Grosbeil ( I 2 OIOI) allowed 
to vary the relative contribution of different groups of orbits 
(binned in energy) in the solution, similar to the Schwarzschild 
method (see below). An advantage of this method is that it uses the 
knowledge of the orbital structure of non-axisymmetric galaxies; 
nevertheless, it still is limited by the assumptions of the functional 
form of potential and perturbation, and has been only used in two 
dimensions. 

A completely different and very generic alternative is of¬ 
fered by combining W-body simulations (which by definition are 
dynamically self-consistent, even though not necessarily station¬ 
ary) with some adjustment method to bring the model towards 
the required properties. The M a de-to- Measure (M2M) method, 
introduced by ISver & Tremaine! ( Il996ll . deals with a fixed, pre¬ 
determined potential, and adjusts the weights of individual particles 
in the course of simulation to match the prescr ibed observable prop¬ 
erties. [Bissantz, Debattista & Gerhard ( 12004) and ide Lorenzi et ^ 
l l2007h extended this method to rotating potentials and included a 
likelihood-based treatment of kinematic observational constraints. 
This implementation also has an option of adjusting the grav¬ 
itational potential together with the evolving density (using a 
spherical-harmonic expansi on), thus makin g the resulting model 
dynamically self-consistent. IPehnenl ( 120091) improved the method 
by introducing a new scheme for adjusting pa rticle weights , which 
has later been incorporated into other codes. iLong & Mad ( I 2 OI 0 I) 
and iHunt & Kawatal ( 1201 3l) developed two more impl ementations 
of M2M method, applicable also to disc galaxies dLong et al.l 
I 2 OI 3 I ; iHunt. Kawata & Martell l2013h . Another class of methods 
based on “guided” A-body simulations are the iterative methods of 
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iRodi onov. Athanassoula & Sotniko'v^ l l2009h and lYurin & Springeil 
l20l4) . In these approaches, the evolution of a live A'^-body system 
is followed for a certain time (an episode), and then the velocities 
of particles are adjusted according to some scheme, to match the 
required properties of the system. The whole process is repeated 
iteratively until a convergence is reached. Both M2M and iterative 
methods are very flexible a nd capable of dealing wi th any geome¬ 
try (although the method of lYurin & SDring^l2014l is restricted to 
axisymmetry); they are however rather costly in terms of computa¬ 
tional demand. 

Finally, the or bit superposition method, introduced by 
ISchwarzschild l ll979h . has been widely used in the last decade to 
create stellar-dynamical models of external galaxies. There exist 


sever al i mplementations for axisy mm etric geometry ([Crettgn eLal 
I 999 I: iGebhardt et alj I 2 OOOI : IValluri. Merritt & Emsellem 


2004) and one for non-r otating triaxial potentials with cores 


I van den Bosch et^l2008h . which are capable of dealing with a 
variety of observational constraints and perform a model search to 
find the best-fitting parameters of the potential and their uncertain¬ 
ties. There are also implementations targeted fo r a more theoretical 
usage ( without observ a tional constra i nts), e . g. iMerritt & FridmanI 
l ll996h : ISiopisI ( Il998l) : iThakur et alj i l2007h : IVasilievI <20131) . All 
these studies considered galactic models that are not too far 
from spherical, i.e. they could be axisymmetric or triaxial, but 
not strongly flattened. Applica tions to barred disc galax i es are 
rather scarce: IPfennige^ (Il984bl) and IWozniak & Pfennigej i ll9971) 
considered two-dimensional models for the motion in the galactic 
plane of a barred potential, and several s tudies jZhaol Il996bl : 
iHafner et alJbOOOl : IWang et al.l[20ll |2013l) concentrated on the 
models of Milky Way bulge, without attempting to construct a 
model for the entire disc. 

In this paper, we present a publicly availablj^ implementation 
of the Schwarzschild method suitable to deal with multicomponent, 
arbitrarily flattened, non-axisymmetric potentials with figure rota¬ 
tion. It is a further development of the computer code SMILE, de¬ 
scribed in detail in IVasili^ (l2013l) . In the present version, it re¬ 
mains “a theorist’s tool”, in the sense that it does not include the 
possibility to deal with observational kinematic constraints; how¬ 
ever, it offers a very generic framework suitable for all types of 
galaxies, and in the future we plan an observational extension of 
the code. 


3 TECHNICAL DETAILS OF SCHWARZSCHILD 
MODELLING 

To construct a self-consistent equilibrium model with the 
Schwarzschild method, one takes the following steps: 

(i) Choose a model for the gravitational potential. 

(ii) Compute a large number (typically > 10'*) of orbits in this 
potential, sampling, in as much as possible, the entire phase space. 
Each orbit is computed for many (~ 100) dynamical times and its 
properties are stored in a discretized way (see Section[33]l. 

(iii) The model is constructed as a weighted superposition of 
these structural blocks (orbits), with the weights being computed 
as a solution to some optimization problem. 

If the method is used in the observational context, then a series of 
models with varied parameters of the potential is constructed, and 
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the likelihood of each model in fitting the observational data is used 
to derive the confidence ranges of the potential parameters. 

We describe these steps in detail in the following sections, 
focusing on the particular features that are necessary to deal effi¬ 
ciently with models of disc galaxies. 


3.1 Potential approximation 


At the very least, a flexible and accurate method for representing 
an arbitrary potential of a highly flattened system is required. Of 
course, one may use a combination of analytical potential-density 
models which approximate the target system, such as a Miyamoto- 
Nagai disc or a Ferrer s triaxial bar, which have been used in man y 
previous studies (e.g.. lAthanassoula et al.|[l983l : IPfennigeill 19843) . 
This approach may well capture important properties of the sys¬ 
tem, but it is difficult to control the systematic errors from para¬ 
metric models. On the other hand, approximating the potential by 
a suitably chosen basis-set expansion is a more general approach, 
although it is not trivial to do it efficiently for a strongly flat¬ 
tened system. Most work that has been done in this direction is 
based on the spherical-harmonic expansion of both density and 
potential up to a given order (max in the angles, and express¬ 
ing the radial de pendence of expansion coefficients as a sum over 
basis functio ns Jciutton-Brock| ll97^ : iHernguist & Ostrik^ll992l : 
IZhadll996al: IWeinberd IT99^. or as explicit fu nctions of radius 
( lMcGlvnijll984 ISellwoodl2003l: IVasili^l2013l) . This technique, 
especially in its last variant, offers a very good approximation 
of the potential of systems which are not too far from spheri¬ 
cal (i.e. for elliptical galaxies), but becomes increasingly ineffi¬ 
cient for disky systems as the order of angular expansion, nec- 
essary to resolve the thinness of the disc, inc reases (for instance, 
[Hollev-Bockelmann, Weinberg & Kat3 (l2005h used (max = 36). 

On the other hand, for strongly flattened systems it is more 
natural to work in cylindrical coordinates R, z, 0 rather than spher¬ 
ical r, 6, (fi. For two-dimensional systems, an efficient basis set 
may be constructed u sing approp riately scaled Bessel functions 
dClutton-Brockl 1 19721 : lOianl Il99^. Th e third dimension may be 
added in several different ways. EarnI Il996h presents a basis se t 
for density models separable in R, z, while iRobiin & Em] ( 1 19961) 
genera lize this approach for other coordinate systems besides cylin¬ 
drical. IWeinbera d 19991) resorts to a numerical solution of the 
Sturm-Liouville equation to obtain the basis functions in R, while 
using harmonic expansion in rj) and 2 ; (the latter in a finite-size 
sheet). [Brown & Papaloizoul d 19981) employ a similar approach 
with both R and 2 basis functions being expressed by suitably 
transformed polynomials, tailored to a given lowest-order density 
distribution. In a special case of a separable axisymme t ric de nsity 
profile in cylindrical coo rdinates, iKuiiken & Dubin^ d 19951) and 
IPehnen & Binnevldl998l) used an approximation based on the split¬ 
ting of the potential into a separable exactly representable flattened 
part and a spherical-harmonic expansion for the remaining weakly 
non-spherical part. Einally, in A^-body simulations often a three- 
dimensional cylindrical g rid is used (e.g. lPfenniger & Friedlill9^ : 
ISellwood & Valluril 19971) : however, it usually employs a low-order 
(linear) interpolation for force computation, which might be insuf¬ 
ficient for high-accuracy orbit integration. 

Given a good experience of working with a non-parametric 
representation of spherical-harmonic expansion coefficients as 
spline functions in scaled radius, we decided to use a similar ap¬ 
proach for flattened systems in cylindrical coordinates. Namely, we 
represent the potential as a sum of Eourier terms in the azimuthal 
angle cj>, with the coefficients of expansion being smooth functions 
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in the meridional plane R, z, that is, two-dimensional cubic splines 
in suitably scaled coordinates. We refer to the appendix for a more 
technical description and accuracy tests of this potential approx¬ 
imation. Overall, its performance is good enough to be used for 
approximating an arbitrary flattened non-axisymmetric mass distri¬ 
bution. Moreover, in multicomponent models we may choose the 
most suitable potential representation for each component, for in¬ 
stance, the cylindrical spline for the disc galaxy and the spherical- 
harmonic spline expansion for the halo. 


3.2 Generation of initial conditions for the orbit library 


The core feature of the Schwarzschild method is that it automati¬ 
cally picks up the orbits that are most suitable for representing the 
target density profile. Nevertheless, the efficiency of the method 
and the very existence of the solution do depend on the procedure 
employed to sample the initial conditions of the orbit library. As 
an extreme example, suppose we wish to create a spherical system 
in dynamical equilibrium using purely radial orbits (even though it 
is guaranteed to undergo a violent radial-orbit instability, this does 
not invalidate the concept of equilibrium solution). If we pick up 
initial velocities with random orientations, for instance using the 
isotropic distribution function, the chances to find a purely radial 
orbit are zero. On the other hand, if we use a spherical Jeans equa¬ 
tion with the velocity anisotropy coefficient /3 = 1 — 
equal to unity (where and at are the velocity dispersion in the 
radial and tangential directions, respectively), then all our orbits 
will be radial. We could also start all orbits with zero velocities, 
and do not even need to seed the initial positions of the orbits to 
follow the required density profile - the solution obtained with the 
Schwarzschild method will converge to it, as long as we ensure 
that there are sufficiently many (not necessarily all) radial orbits in 
the orbit library. As another example, using isotropic velocities in 
a rotationally-supported disc galaxy would be much less efficient 
than, for instance, launching orbits with some small velocity dis¬ 
persion about the mean circular velocity; however, if this velocity 
dispersion is too small for a given thickness of the disc, a self- 
consistent model could not be constructed either. 

This suggests that the choice of initial conditions is an im¬ 
portant ingredient of the model, and should, in as much as pos¬ 
sible, correspond to the target system to be constructed. In the 
traditional approach, used in nearly all previous studies (e.g. 
iMerritt & Fridman|ll996l ; Ivan den Bosch et aDl2008h . one takes a 
small number of discrete values of energy E and populates the 
initial conditions fr om several “start-spa ces” at each energy. Sta¬ 
tionary start-space ( l^hwarzschildll 1993h consists of points cover¬ 
ing the equipotential surface with zero initial velocities; principal- 
plane start-spaces contain points with one coordinate (say, z) and 
two complementary velocity components (Vx,Vy) being zero, and 
the remaining component Vz assigned to yield the required total 
energy. In the axisymmetric case, it is sufficient to sample the po¬ 
sitions in the meridional plane on the zero-velocity curve for the 
given values of energy and angular mo mentum (e.g. Crgtton et al.l 
I l999h . For models with figure rotation. ISchwarzschildl (Il982l) pro- 
posed a y — a start space, with initial conditions taken along the 
intermediate (y) axis and with velocity directed at the angle a in 
the X — z plane. 

All these methods place initial conditions on some sort of reg¬ 
ular grid, which may not always be a good choice. For instance, 
with 10^ orbits we only have ~ 20 grid nodes per each dimension 
(the grid of initial conditions for triaxial models is usually three- 
dimensional), which leads to excessive granularity of the orbit li¬ 


brary. As shown bv IVasiliev & Athanassoulj l l2012h . such models 
are not in perfect equilibrium and demonstrate a noticeable evolu¬ 
tion in the course of A'-body simulations at early times. An alterna¬ 
tive approach, used in the aforementioned paper, is to assign initial 
conditions randomly, for instance, sampling them from an isotropic 
distributio n function computed with the Eddington inversion for¬ 
mula (e.g. ISinnev & Tremainell2008L equation 4.46) for a spheri¬ 
cal potential-density pair that approximates the real, non-spherical 
model. In this approach, the positions of particles are first seeded 
in accordance with the true three-dimensional density profile of the 
model, then the velocities are assigned by drawing them from the 
Eddington distribution function. 

In the present implementation, we introduced several addi¬ 
tional methods for generating the initial conditions for velocities, 
aside from the Eddington sampler. The first two are based on the 
Jeans equations for the spherical and axisymmetric systems. The 
former is rather trivial, and its only improvement over the Edding¬ 
ton sampler is the facilitation of creating models with strong radial 
or tangential velocity anisotropjQ. The axisym metric Jeans equa - 
tions are used in the anisotropic formulation of lCappellaril ( l2008l) : 
the velocity dispersion ellipsoid is assumed to be aligned with the 
cylindrical coordinate system, and the anisotropy coefficient in the 
meridional plane /Jm = 1 — is assumed to be constant. 

These assumptions are rather strong and moreover they are not sat¬ 
isfied for realistic galaxies; for an alternative method of solving 
the axisymmetric Jeans equati ons under assumpt i on of radial align¬ 
ment of velocity ellipsoid see lYurin & Sprint ( l2014ll . Neverthe¬ 
less, for the purpose of generation of initial conditions this simple 
choice is sufficient. Aside from the anisotropy coefficient /Im, the 
other free parameter is the degree of rotation support k\ the mean 

_ /- -N 1/2 

azimuthal streaming velocity is given hy = k ^r) 

With this definition. A: = 0 corresponds to no net rotation and 
k = 1 corresponds to equal velocity dispersions in R and 0 di¬ 
rections. 

Another approach, which could more tailored to the 
particular properties of the potential, is ba sed on the knowl¬ 
edge of periodic orbit fam il ies ( e.g. _ Athanassoula et all 

19831; Con topoulos & Gro sbdil 19891; Athanasso ulal 199? 
Hasa n, Pfenniger & Normanll 19931 ; Skokos, Patsis & Athanassoulal 


I 2 OO 2 I 1 . First. one constructs a “periodic orbit map”, which consists 

of the initial conditions for various periodic orbit families (most 
importantly, prograde loop orbits) in the range of values of Jacobi 
constant Ej from the bottom of the potential well to the corotation. 
Then one seeds initial conditions for the actual orbit library by 
picking them from the map and adding a random perturbation to 
the velocity components, to diversify the orbit library. We have 
experimented with this approach, but did not find a noticeable 
improvement over the simpler method based on the axisymmetric 
Jeans equations; therefore, in the rest of the paper we use only the 
latter. 


^ We do not argue that the Jeans equation approach is supe¬ 
rior to the distribution-function-based method, and it is known that 
the former may lead to significant biases in some applications 
I Kazantzidis. Magorrian & MoorejEoo4) . However, in the Schwarzschild 
method the distribution function is recovered numerically, and thus these 
drawbacks are irrelevant. 


© 2015 RAS, MNRAS 450, 2842-2856 














































2846 E. Vasiliev, E. Athanassoula 


3.3 Obtaining the self-consistent solution 

After computing the entire orbit library, the next step in the 
Schwarzschild method is to obtain the weights of orbits, such that 
the weighted superposition of density distributions of each orbit re¬ 
produces the target density profile of the model. Additionally, one 
often requires that some other constraints be satisfied, for instance, 
fitting the line-of-sight velocity profiles from the model to the ob¬ 
served values. The present implementation allows only for very 
rudimentary kind of kinematic constraints: one can specify the pro¬ 
file of velocity anisotropy coefficient as a function of radius. We 
will not use this feature for the models of discs, and will set /3 = 0 
for the haloes. 

The density profile of the model can be represented in a dis¬ 
cretized way either as an ar r ay of masses of cells in a spatial grid, 
or, as suggested in lVasilie^ (l2013h . as an array of coefficients de¬ 
scribing the basis-set or spline spherical-harmonic expansion of the 
potential corresponding to this density. Since spherical-harmonic 
expansions are poorly suited for highly flattened systems, we use 
the former, traditional approach of partitioning the configuration 
space into a three-dimensional grid, computing masses of each grid 
cell rUc by integrating the target density profile over the volume 
of the cell, and recording the fraction of time toe that each orbit 
spends in each cell. For elliptical galaxies, a grid based on concen¬ 
tric shells in radius and a certain partitioning in angles is typi cally 
used fe.g. lMerritt & Fridmanll996l : lvan den Bosch et alj2008h . For 
flattened systems we instead use a grid aligned with cylindrical co¬ 
ordinates (covering the region 0 < R < i?max,0 < a < Zmux 
with Nr X Nz rectangular cells, and split into cells in the az¬ 
imuthal direction, equidistantly in (j)). The grid nodes in the merid¬ 
ional plane are assigned in such a way that the total mass in each 
slab dR d(j> dz Rp{R, (j), z) and a simi¬ 

lar expression for ) is the same for each slice in z direction or 
cylindrical shell in R direction. Of course, if the density profile is 
not separable in R, z then the masses of each cell are not equal, but 
they usually vary only by a factor of few. 

The self-consistent solution for orbit weights Wo is given by 
the matrix equation Wotoc = rrio, c = 1... A^ceii- This is 
not just a standard linear-algebra problem, as we require the so¬ 
lution vector Wo to be non-negative; in this formulation it is called 
a linear programming problem. In practice, an exact solution may 
not always be feasible; moreover, it is worthwhile to apply some 
sort of non-linear regularization to improve the smoothness of the 
solution. A suitable formulation that meets these requirements is 
the quadratic programming problem: minimize R = Rceii + Rorb 
while keeping Wotoc — rrio = So- Here So is the deviation in 

the mass of c-th cell from its required value. The penalty function 
R is split into two parts: one is responsible for minimizing the devi¬ 
ations from the self-consistent solution (which would have So = 0), 
the other regularizes the orbit weights and optionally imposes addi¬ 
tional penalties for using particular types of orbits (e.g. orbits with 
negative 2 -component of angular momentum, or chaotic orbits). 

For the first task, several methods have been suggested in the 
past (here a... are adjustable parameters): 

(i) J^constr = ai I] |(5c| (e.g. |Siopis|| 19981) ; 

(ii) J^constr = Q2 y~) (least-SQ uare minimization, used in 
iMerritt & Fridmadl 199d ; lzhaoll 1996bl) ; 

(iii) J^constr = I]{0if|5c| < otorUo, otherwise 00 } (con- 
st raints should be satisfied w ithin a given fractional error ao, as 
in ivan den Bosch et al.ll2008h . 

The first method is easily implemented in the context of linear op¬ 


timization problem by introducing 2A^ceii additional non-negative 
variables po, Vo such that So = po — Vo', thus the penalty function is 
just proportional to the sum of these variables (of which at most one 
is non-zero for each constraint). Similarly, in the second method, 
which traditiona lly has been solved using the non-negative least- 
squares method l lLawson & Hansoi]|l974h . the introduction of the 
same additional variables transforms it to a quadratic optimization 
problem with J^constr = 0.2 X]c(Fc + which we have found 
to be more numerically efficient. Finally, the last method can also 
be reformulated as a linear optimization problem with twice larger 
number of equations. All these methods are available in SMILE; we 
usually use a combination of first two methods, so that small devia¬ 
tions in cell masses are penalized linearly and large - quadratically. 

The regularization method that we use for the second task just 
aims to achieve an uniform distribution of orbit weights: 


Rorh — 


. ^orb 

0=1 


where Wo is the weight prior (in the case of uniformly sampled ini¬ 
tial conditions it is just the mean weight of an orbit in the model), 
and A is the regularization coefficient. Adjusting the values of ai ^2 
and A, one may vary the relative severity of constraint violation 
versus irregularities in the distribution of orbit weights. In practice, 
for a well-behaved model with a sufficiently large number of orbits 
(Norb S> A^ceii). we expect all constraints to be satisfied exactly, 
thus only the second term in R remains nonzero. However, as we 
have discovered while working with potential approximations con¬ 
structed from A^-body snapshots, sometimes the mass of several 
cells could be computed from the density model with a rather large 
error and cannot be satisfied by the Schwarzschild solution. In this 
case it may be better to sacrifice a few cells in order to keep the 
model reasonably smooth. 

An extension over the previous version of SMILE is the full 
support for multicomponent models. This is implemented as fol¬ 
lows: the orbits are computed in a combination of several po¬ 
tential models (e.g. triaxial Ferrers bar plus an exponential disc 
plus a spherical Hemquist halo), and the solution of the self- 
consistent problem uses an arbitrary subset of these density mod¬ 
els (e.g. the first two for the “disc” component and the last for 
the “halo” component, which are constructed independently and 
may use different methods for seeding the initial conditions). 
Whil e multicomponent Schwarzsch ild models have been used be¬ 
fore dCapuzzo-Dolcetta et al.ll200% . the present implementation is 
much more general and versatile. This extension required some 
modifications in the Jeans and Eddington initial condition samplers, 
to distinguish between the total potential of the entire model and the 
intrinsic density profile of a particular component. 


4 TESTS 

4.1 Axisymmetric disc-bulge-halo models 


In the first test, we construct composite models of an ax¬ 
isymmetric disc galaxy with a bulg e and a halo, using several 
differ ent methods: MaGalie jBoilv. Kroupa & Penarrubia 
l200lh . a desce ndant of Hernquist’s (1 993) BUILDGAL code; 
MKGALAXY I M^illaii & DehnenI l2007h : GalactICs 

dWidrow, Pym & Dubinski||2008 h, whic h is the la test version of the 
metho d of iKuiiken & Dubinskil ( Il99^ ; GALIC dYurin & Springell 
l2014h ; and SMILE. All these c odes are public ly available: the first 
two as part of NEMO toolbox dTeuberj[l995l) . the third as part of 
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Figure 1. Circular velocity curve of the axisymmetric disc-bulge-halo 
model used in the test. 

AMUSE framework dPelupessv et aklEoiSh . and the last two at the 
authors’ websites. 

The first complication that we encountered is the different 
specifications of galaxy parameters in these codes. All methods 
employ various parametrization of three components - disc, bulge 
and halo density profiles, with fixed but sometimes incompatible 
functional forms, and in several cases, specified only indirectly and 
non-intuitively: for instance, instead of component masses GALAC- 
TICS uses the depth of gravitational potential, and the disc scale- 
length in GALIC is determined by its mass and an obscure halo 
spin factor. After some experimentation, we have selected the fol¬ 
lowing parameters, which produced nearly identical mass models 
for all codes: an exponential disc with density profile paisc = 
Mdisc/(47ri?oZo) exp(—7?/i?o) sech? {z/ zq)-, a bulge with scale 
radius ruuige and either Sersic profile with the index n = 2, or 
Hernquist profile; and a Navarro-Frenk-White (NFW) halo with 
cutoff at 10 scalelengths, or an equivalent Hernquist halo. The 
masses and sizes of these components are: Maisc = 1, iZo = 1, 
zq = 1/8, Mbuige = 0.25 (for the Sersic profile), ruuige = 0.2, 
AThaio ~ 25 (depending on the functional form of cutoff), and halo 
scale radius of either 5 (for NFW) or 10 (for Hernquist). The ro¬ 
tation curve of this model is shown on Fig. [T] The choice of disc 
thickness implies a moderately warm disc, and the masses and scale 
radii of spheroidal components give them a rather large contribu¬ 
tion to the total gravity: both factors aim at avoiding disc instabil¬ 
ities, but at the same time stress the necessity to take into account 
the gravity of all mass components when constructing the disc dis¬ 
tribution function. 

Exponential disc is the standard option for all these codes, but 
the density profile of the spheroidal components is constructed dif¬ 
ferently: MKGALAXY and GALACTICS imply that the density pro¬ 
file follows the equipotential surfaces, which are somewhat flat¬ 
tened in the central parts. Consequently, the zjx axis ratio of the 
bulge is close to 0.8 in these codes, so we adopted this value for the 
other methods. The halo axis ratio varies from ~ 0.7 in the centre 
to about 0.8 at 7Z = 2 (the maximum of disc rotation curve) and 
then approaches unity as R increases further. For SMILE we have 
incorporated this spatially variable flattening into the input density 
profile; the other two codes do not allow variable flattening and we 
have assumed a spherical halo for them. The velocity dispersion of 
bulge and halo was required to remain isotropic. Most differences 
between codes are manifested in the velocity distribution of the disc 
component. 

The models constructed with these five methods all had 
(160, 40, 800) X 10® particles in the disc/bulge/halo components, 
and were checked for stability by evolving them for 100 time units 



0 1 2 3 4 5 

R 

Figure 2. Velocity dispersion profiles of the disc component, for five 
different methods. Shown are values of an (dashed line), a^j, (dotted 
line) and cr^ (dot-dashed line), normalized by the value of veifical veloc¬ 
ity dispersion for an isolated isothermal exponential disc: o'z,iso(l?) = 

Mdiac zo/{2Rl) exp(-i?/(2Ro)), where M^iac = 1, fio = 
1, 20 = 1/8 are the disc mass, scale radius and scaleheight. 

Thicker lines are the initial models, thinner lines with shorter dashes ai‘e the 
models evolved for 100 time units. 

Also shown (in solid lines) are radial velocity profiles corresponding to the 
value of Toomre parameter Q = 2 (in the second panel) and to the solution 
of axisymmeti'ic Jeans equation with /?m = 0.6 (in the last panel, served as 
initial conditions for orbit library). 
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Figure 3. Toorare Q parameter for the disc component of the axisymmetric 
models built with different methods. 


with the A'^-body code GYRFALCON (lDehnenl20Q(ll20Q3) : the soft¬ 
ening length was set to 0.01 for disc and bulge and 0.03 for halo 
particles, and the timestep was 2“^. Fig. shows the radial pro¬ 
files of three diagonal components of velocity dispersion tensor in 
cylindrical coordinates - thicker lines for the initially constructed 
models,and thinner lines for the evolved models. The velocities are 
rescaled in such a way that the vertical velocity dispersion of an 
isolated isothermal exponential disc would appear as a constant line 
equal to unity. 

First thing to note is that the codes produce a diversity of ve¬ 
locity profiles, determined by different constraints. To assign the 
radial velocity dispersion, MAGALIE and MKGALAXY use Toomre 
Q parameter, GALACTICS assumes its exponential dependence on 
radius, and GALIC and SMILE solve axisymmetric Jeans equation. 
The vertical velocity dispersion is simply taken from an isolated 
isothermal disc in MAGALIE and MKGALAXY, resulting in a too 
low value in the centre and in the outermost part of the disc (where 
there is a substantial contribution to the gravity from the bulge and 
the halo, correspondingly). SMILE and GALIC take this contribution 
into account, while GALACTICS is halfway between these options. 
GALIC enforces — an, while other codes determine from 
other considerations, resulting generally in > (jr. We note that 
the initial conditions for the orbit library in SMILE are assigned with 
= (^R = (Jz/x/l ~ /3m (with the choice of /Jm = 0.6), however 
the actual orbit integration results in a different velocity structure, 
which we do not constrain. Adding such constraints is not difficult, 
but as Figs [3 and (solid line in the bottom panel) show, neither 
a constant Q nor a constant (jr/cjz are good choices for the tar¬ 
get velocity dispersion profile, thus raising the question of the most 
natural requirements for it. 

The results of the simulations are also quite diverse. The mod¬ 
els that assumed the vertical velocity dispersion to be equal to that 
of a self-gravitating isothermal disc (MAGALIE and MKGALAXY) 
have quickly increased it in the central parts, while the disc thick¬ 
ness correspondingly dropped in the centre. In all models, the ver¬ 
tical thickness and velocity dispersion of the disc component has 
slowly but steadil y grown over tim e, presumably as a result of two- 
body relaxation l lSellwooal2013l) . and due to slowly developing 
disc instabilities. MAGALIE and GALIC did not preserve the ini¬ 
tial profiles of Or and (the former code - only in the outer 
parts). The biggest limitation of the latter code stems from enforc¬ 
ing gr = (T0 (or, in other words, assuming the streaming parame¬ 
ter k = 1): while there are other possible choices in the code, none 
corresponds to the actual variation of k seen in other models. It 
might be better not to fix the value of in this approach. GALAC¬ 
TICS generated a model reasonably close to equilibrium, and SMILE 


fared best in this aspect. Note that in the latter model the velocity 
dispersion profile is considerably higher at large radii than in other 
methods. In terms of computer resources, MAGALIE and GALAC¬ 
TICS take a few minutes to create the models, MKGALAXY takes 
about an hour, SMILE takes a few CPU hours (for the total number 
of 5 X 10“^ orbits in our models). GALIC took about 100 CPU hours 
for 100 iteration steps, but since the model taken after 10 iterations 
is essentially the same as after 100, a more truthful execution time 
is about 10 CPU hours. GALIC used ~ 32 Gb of memory, while 
SMILE needed about 1 Gb, and the other codes much less. 

The results of this test demonstrate that SMILE is capable of 
creating multicomponent systems in almost perfect equilibrium; a 
big improvement over previously used codes. It is important to note 
that, unlike other codes considered here, it allows an arbitrary den¬ 
sity profile for (any number of) mass components, although in the 
present version it does not impose any constraints on the velocity 
structure. 


4.2 Rotating triaxial Dehnen model 


As a second test, we attempt to construct models of mildly flat¬ 
tened, cuspy density profiles, typical for elliptical galaxie s, having 
a vary ing degree of figure rotation. We consider a 7 = 1 IPehnerJ 
( I1993I) model with axis ratios x : y : z = 1 : 0.75 : 0.5, which ro- 
ta tes about its short axis. Th i s is on e of the models that was studied 
in lPeibel, Valluri & Merritj ( 1201 ih in the context of orbit analysis. 
The distinct feature of moderately cuspy (7 = 1) triaxial Pehnen 
models is the rich network of resonant and thirFI orbits, which are 
clearly s een on frequency map plots (e. g. IValluri & Merritllll998L 
fig. 9, or IVasiliev & Athanassoulall2012L fig. 2). Orbits associated 
with these resonances have a variety of shapes and are important 
ingredients of self-consistent models, replacing regular box orbits 
that appear in models with constant-density cores (the non-resonant 
box orbits are generally chaotic in cuspy potentials). The effect of 
rotation o n regular box and thin orb its is the so-called “envelope 
doubling” dde Zeeuw & Merritt! 1983ll . caused by changing the sign 
of the Coriolis force as the orbit travels with positiv e or negative 
angular momentum. IPeibeL Valluri & Merrid l l201lh have shown 
that this effect does not stabilize the chaotic orbits by converting 
them into regular o nes that avoid passing thro ugh the centre, as has 
been suggested bv iGerhard & BinnevI l ll985h : on the contrary, the 
fraction of chaoti c orbits in rota ting models is higher than in non¬ 
rotating ones tcf. lMuz^l2006h . and increases with pattern speed, 
although some of them become converted to regular tube-like or- 
bit s when the rotation is high eno ugh. Thus, the question raised 
by IPeibeL Valluri & Merritl ( 1201 ih is whether the self-consistent 
models could be constructed for moderate to high pattern speeds, 
and we address this question below. 

We take several values for the pattern speed (see table 1 in that 
paper), corresponding to the corotation radius Rc equal to 20, 10, 
5 and 2, and construct Schwarzschild models with 5 x 10"^ orbits, 
~ 10® constraints, and the initial conditions generated using ei¬ 
ther Eddington or axisymmetric Jeans equation (the latter produces 
more orbits that rotate in the same sense as the figure of poten¬ 
tial, but the results did not differ significantly between these two 


® A thin orbit is confined to a two-dimensional, possibly self-intersecting 
sheet in the configuration space; its fundamental frequencies uji of motion 
in three directions satisfy a resonant relation X]i=i = 0 with integer 
ai. 
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Figure 4. Axis ratios of four rotating Dehnen models considered in Sec- 
tion l4.2l labelled by their corotation radius. Thick lines are the initial mod¬ 
els, and thin lines show models evolved for 1000 time units. The original 
models had a variable y : x axis ratio, ranging from 0.75 in the central part 
to 1 beyond the corotation radius. Evolved models largely maintained their 
initial shape, although becoming somewhat rounder (this effect was most 
noticeable for models with extreme angular momentum). 



Figure 5. Evolution of pattern speed of four rotating Dehnen models, nor¬ 
malized to its initial value (which equals 0.0109, 0.0287, 0.0745, 0.257 for 
models with corotation radii Rc = 20,10, 5, 2, correspondingly). 


approaches for these mildly flattened models). It was quickly dis¬ 
covered that the model cannot maintain elongated shape beyond 
corotation, so we modified the density profile to have a variable 
X : y axis ratio, which was kept at the value 0.75 up to 0.67?c, 
and further out gradually changed towards unity (Fig.|4j; further¬ 
more, we did not attempt to constrain the azimuthal distribution of 
mass outside i?c- We confirmed that the degree of chaoticity of all 
major orbit families in general increased with the pattern speed. 
Considering the orbit population in the region R < 2, which en¬ 
closes roughly half of total mass, we found that all models except 
the most rapidly rotating one contain around 40% short-axis tubes 
and 20% resonant and thin orbits, while in the most rapidly rotating 
model nearly 80% orbits are short-axis tubes. 

We then evolved W-body representations of each model for 
1000 time units with GYRFALCON, using 5 x 10® particles and a 
softening length of 0.01. We measured the shape and orientation of 
the models, using both Fourier analysis of the sur face density and 
the method of inertia tensor (see lZemp et alj201 iL for a discussion 
of its variants). The models turned out to be rather stable: by the 
end of the simulation all of them largely retained their shape and 
pattern speed (Fig. |3, the latter changed by at most 15% (recall 
that the most rapidly rotating model has performed around 40 full 
revolutions by that time). 

We also briefly explored how the pattern speed is related to 
the angular momentum of our models. With the default settings, 
the total angular momentum turned out to be roughly proportional 
to the pattern speed. For the model with = 5 we created two ad¬ 
ditional solutions that minimized and maximized the total angular 
mo mentum, adding co rresponding terms into the penalty function 
(e.g. |Pfennigen[l984bl) . The former resulted in a model where the 
streaming velocity and pattern speed were directed in the oppo¬ 
site sense, while in the latter case almost all orbits rotated in the 
same direction (Fig. [Hi. Both of these models were confirmed to 
be reasonably stable in W-body simulations; the model with maxi¬ 
mum streaming motion actually had increased the pattern speed by 
~ 15%. The axis ratios of these models by the end of simulation 
were closer to unity than for the “default” models, although they 
still retained distinctly triaxial shapes. 



Figure 6. The rotation parameter Ar jEmsellem et alj|2007h for triaxial 
Dehnen models with figure rotation, seen along the intermediate (y) axis, 
plotted as a function of distance to the centre r = + z'^', half of 

the mass lies within r ^ 1.4. This parameter characterizes the amount of 
streaming motion, normalized to the total mean-square velocity; a value of 
0.1 separates slow from fast rotators. In the “defaulf’ models with no kine¬ 
matic constraints, streaming motion is generally proportional to the pattern 
speed, but there is a lot of freedom in it: for two additional models with the 
same pattern speed but maximizing or minimizing the angular momentum, 
it varied in a much wider range (the latter actually was counter-rotating with 
the figure of potential, which is shown by a negative value of Aj^). 


4.3 Barred disc galaxy model 


For the next test, we consider an analytic potential model with a 
Miyamoto-Nagai disc, a Ferrers bar and a Plummer bulge, with 
the s ame parameters as in ISkokos, Patsis & Athanassoulal 

(|2002|); sim ilar m ode ls were previously exp lored _ in 

Pfenn igeJ (Il984ah: [Hasan, Pfen niger & Norman (Il993l): 
Patsis. Skokos & Athanassoulal 320021) : iManos & Athanassoula 


120111 etc., in the context of orbit analysis, while jpfl 


ennigei 


1984bl) has constructed two-dimensional (in the equatorial plane) 
Schwarzschild models for a similar density profile (without 
the bulge component). The pattern speed is chosen so that the 
corotation radius is just heyond the end of the bar. We combined 
the density profiles of all three components into one model and 
did not require the orbits to fit each of them separately, as in 
Section 14.11 We considered several choices for the number of 
orbits (up to 10®) and density constraints (up to 2 x 10®), as well 
as for the generation of initial conditions (using axisymmetric 
Jeans equation with different choices of Pm, and launching orbits 
from the vicinity of stable periodic orbits). 

The results of all these experiments suggest that there is some¬ 
thing unrealistic about the density model that we attempt to con¬ 
struct self-consistently. While the models with a coarse spatial grid 
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Figure 7. Distribution of orbit weights in a Schwarzschild model of a com¬ 
posite model from Section 1431 as a function of orbit index (sorted in en¬ 
ergy). In the ideal case, it should be narrowly spread around 1/Na^\y = 
10~®, but in this case this happens only at small or large energies; around 
the energy corresponding to the end of the bar the model barely manages 
to satisfy the constraints by using only a few orbits with relatively high 
weights. This usually indicates troubles in finding a self-consistent solution. 


could be constructed, the distribution of orbit weights in the solu¬ 
tion turns out to be very uneven, especially in the region around 
the end of the bar (Fig.|7](. Furthermore, a detailed inspection of the 
surface density map of the Schwarzschild model revealed that it did 
not exactly match the analytic model, having sharper edges at the 
end of the bar, even though their discretized versions were equal. 
An attempt to refine the spatial grid resulted in the model becoming 
infeasible, no matter how many orhits we tried. On the other hand, 
when we restricted the orbits to a two-dim ensional plane and fitted 
the surface density, as in IPfennigetl l ll984hl) . the model was feasible 
and matched the projected density very well; however, with zero 
velocity in 2 direction it is not in a dynamical equilibrium. 

The model also turned out to be unstable when tested by an 
A^-body simulation and quickly (within one rotation period) trans¬ 
formed itself into a rounder and more slowly-rotating shape, with a 
less prominent bar. We repeated the experiment with a model that 
lacked a central bulge, and it also demonstrated a similar slowdown 
of pattern speed (by some 30% after five rotation periods), while the 
bar became shorter and narrower. We thus rule out the possibility 
that the bulge was the reason for this evolution. 


4.4 Reconstructing W-body snapshots of barred disc models 
from their density profiles 

In the previous sections, we used analytically defined potential- 
density pairs to create our Schwarzschild models, while in this sec¬ 
tion we will use densities and potentials extracted from A-body 
simulations. The latter approach is, of course, much more realistic 
than the former, because the components are fully self-consistent 
instead of rigid, and, furthermore, because they have grown natu¬ 
rally under the influence of the total gravity of their system. Thus 
their shapes and mass distributions can be much more complex than 
what simple analytic formulae can describe. 

As simulations have sho wn, barred galaxies are not steady- 
state systems (for a review see lAthanassoulaluOlSl and references 
therein). The bar region is ready to emit angular momentum, 
while the resonant regions in the outer disc and particularly in 
the halo are ready to absor b it iLvnden-Bell & KalnaisI 19721 : 
iTremaine & Weinberd 1 19841 : lAth anassou iT^ lTOOl | 20 ^ so 


that there is a redistribution of angular momentum within the 
galaxy, as a re sult of which a number o f the properties of the 
galaxy change jLittle & Carlbergiri991alf3 : lOebattista & Sellwood 


2000|^ Athanassoula&Misirioty_ 25oJ_ Athanassoulal l2003l : 


O’Neill & Dubinski 20031 : 1 Valenzuela & Klvpinl 2003 ). Thus bars 
become longer and stronger and their pattern speed decreases with 
time. 


Thus barred galaxies will necessarily have some non- 
stationarity and this goes against the assumption of a steady state (at 
least in the rotating frame), implied by the Schwarzschild method. 
They therefore present an additional challenge for this method. Can 
such models be built? And, if yes, can their evolution match that of 
the density distribution we are trying to model, or at least provide 
some useful information about secular evolution? These interesting 
questions will the subject of a future paper, so that only a short sum¬ 
mary of the results relevant to the Schwarzschild modelling will be 
given here. 


Although it is now well established that barred galaxies 
evolve, it is still unclear how strong this evolution is. Simulations 
have shown that the extent to which the angular momentum is redis¬ 
tributed, and therefore how important the corresponding changes in 
the galactic properties are, depends on the properties of the simu¬ 
lated galaxy. It is thus difficult, if at all possible, to use simulations 
to set constraints on evolution. Furthermore, observations show no 
consensus, si nce their results advocate evolution rang ing from rel¬ 
atively little (Perez, A.guem_^Mendez-Abreu| 201_^ to quite im¬ 
portant (e.g. iKorm endv & Kennicuttll2004l:l Elmegreen et al.ll2007l : 


ISheth et all2008l:ICheung et al. ll2013l:|Kim et alj20ll) . It is there¬ 

fore necessary to consider both weak and strong secular evolution 
models. 


We will here attempt to model two examples of barred galax¬ 
ies, extracted from simulations chosen so as to have as different 
secular evolutions as possible; one being very little (model A) and 
the other quite strong (model B). These two simulations are similar 
to the MD and MH ones considered by I Athanassoula & MisiriotisI 
(l 2002 h . having the same initial mass distributions for the halo 
and disc, i.e. the same rotation curves (see fig. 1 of that paper) 
and the same disc-to-halo mass ratio. Thus the halo of model B 
is much more centrally concentrated than that of model A. We 
take a snapshot from each of the two simulation at a time very 
near the beginning of its secular evolution phase and apply the 
Schwarzschild method to construct a steady-state, two-component 
dynamical model, using as constraints only the density distribution 
of both the disc and the halo. These are derived in a non-parametric 
way from the particle positions in the original snapshot. For the disc 
we also need to include figure rotation, i.e. we need an estimate for 
the pattern speed. Since we want to test whether it is possible to 
build models in the absence of kinem atic information, we cannot 
use the ITremaine & Weinberd ( Il984l) method to obtain the pattern 
speed rib- Instead, we tried values in a range encompassing the 
“correct” one, i.e. the value obtained from the simulations. 

The original W-body models had Adisc = 2 X 10® disc 
particles and Ahaio = 10® halo particles. In order to reduce 
the no ise and transi e nt fea tures in the disc component, we fol- 
lowed lAthanassoulal ( l2005h and llannuzzi & Athanassoulal ( 1201 5l) 
and stacked five consecutive snapshots from the original simula¬ 
tion, after rotating them so as to align their bar major axes. The 
potential of the disc was represented with a cylindrical spline ex¬ 
pansion with mmax = 6 azimuthal harmonics. For the halo we used 
a spline spherical-harmonic expansion, keeping only the monopole 
term for model A (which was found to be nearly-spherical), and 
additionally the quadrupole ter ms for model B (which is mildly 
triaxial in the centr al parts, see IColfn, Valenzuela & KlypinllTOOfl : 
lAthanassoulall200% . We used 10® orbits for each component, with 
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the disc and halo density model discretized into ~ 1000 and ~ 500 
cells, respectively. 

The library of orbits for the halo component were created us¬ 
ing the isotropic Eddington distribution function computed for a 
spherical approximation of the true density profile, while for the 
disc component we used the axisymmetric Jeans equations, with 
the meridional anisotropy parameter = 0.5. A comparison 
of the latter with the kinematic properties of the disc in the orig¬ 
inal simulation snapshot shows that the assumption of constant /Jm 
is clearly not very good; it varies from /3m — 0.2 in the centre 
to /3m — 0.9 at large radii. Even so, we were able to construct 
Schwarzschild models for both cases A and B. These have a kine¬ 
matic structure which is very similar to that of the original model, 
if the pattern speed is the “correct” one. Apparently, for these mod¬ 
els there is not much freedom in the solution, so that even without 
any additional kinematic constraints (apart from adopting the “cor¬ 
rect” pattern speed value) they turned out to be close to the original 
models. 


The feasibility of constructing a Schwarzschild model does 
not imply anything about its stability properties, or about the 
amount of secular evolution it will have once it is allowed to evolve. 
Note that this angular momentum exchange, and therefore evolu¬ 
tion, can occur only if both the disc and the halo are represented as 
“live” particles, as opposed to a static potential. 


To test the evolution, we converted the Schwarzschild mod¬ 
els to A-body representations using the same number of parti¬ 
cles as the original snapshots. We then evolv ed them using two 
simulation codes: GADGET-2 ( ISpringelll2005ll and GYRFALCON; 
the re sults were quite similar (cf. Fortin. Athanassoula & Lamberj 
[2onh . except that the GADGET snapshots need to be re-centred. 
Below we use the data from the GYRFALCON simulations. The 
simulations were run for 800 time units (same units as in 
I Athanassoula & Misiriotisll2002h : the original models performed 
about 12 full revolutions during this time. 


We first tested the internal stability of each disc compo¬ 
nent by performing simulations following the evolution of each 
disc in the static potential of its host rigid halo, represented by 
a spline-interpolated spherical-harmonic expansion. For this, we 
used the external potential feature in GYRFALCON with the poten¬ 
tial computed using the SMILEPOT library, essentially the same as 
in the original Schwarzschild model. The absence of any slowdown 
(Fig-in dotted lines) confirms that the models are indeed stable. 


The second test involves the full, two-component system. 
We first evolved the simulation snapshots and give the results for 
their pattern speed in Fig. [8] (solid lines). That of model A did 
not substantially change over the course of the simulation, while 
that of model B decreased roughly by a third. We then evolved 
the Schwarzschild models and found a surprisingdiversity of be¬ 
haviour. Models with pattern speeds higher than the “correct” value 
(i.e. the value obtained from the simulation) tend to slow down 
more rapidly, especially in the case A. Models with pattern speeds 
lower than “correct” evolved less, but in opposite directions for the 
two cases. Models with the “correct” pattern speed evolved rather 
little. Thus, the evolution of model A is very similar to that of the 
snapshot it was built to represent. However, we found also little 
evolution for model B, i.e. this Schwarzschild model evolved much 
less th an the snapshot it was built to represent. From lAthanassouQ 
ll2002h . we expect this difference to be linked to the properties of 
the halo rather than the disc. To test this, we created a “hybrid” 
model with the disc taken from the original simulation, and the 
halo created with the Schwarzschild method. This model showed 



0 100 200 300 400 500 600 700 800 

time 


Figure 8. Pattern speed of model A (top) and model B (bottom) as a func¬ 
tion of time. We show the curves corresponding to the original model (solid 
lines), and to Schwarzschild models with different choices for the pattern 
speed (dashed and dot-dashed lines); dotted line shows the model evolved 
in a fixed halo potential. 


little evolution, similar to the purely Schwarzschild model. We are 
examining the reasons for this behaviour in a separate study. 


5 DISCUSSION 

5.1 Advantages and shortcomings of the Schwarzschild 
method 


In this paper, we have used the Schwarzschild orbit superpo¬ 
sition method to construct several quite different self-consistent 
equilibrium galactic models. This method is very general and 
powerful, but one has to keep in mind several considerations. 
First, the basic assumption behind this and virtually every other 
method for constructing dynamical models is that the system 
under study is in a stationary state. This needs not be the 
case, for instance, if we model non-axisymmetric disc galax¬ 
ies, in which a complex interplay between local and global in¬ 
stabilities, resonances, dynamical friction and other mechanisms 
leads to constant evolution, even in the absence of dissipa¬ 
tion. Diffusion of chaotic orbits may also pla y a role in driv¬ 
ing the secular evolution of both elliptical (e.g. Mer ritt & Va lluril 


Il996l : [^ndrup & Siop ijgodj: Vasi liev & Athanassoula||20l A and 

disc galaxies (e.g. Vogli s. Stavropoulos & Kalapotharakoslliood ; 


IContopoulos & Harsoulall2013h . Thus the models that we intend to 
be in a steady state, may turn out to evolve considerably, under¬ 
mining our assumption. In application to real galaxies, it might be 
possible to create a model that matches all existing observations, 
but we cannot ensure that its evolution will be the same as for the 
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original galaxy. The example of Section 1441 shows that we could 
construct a nearly stationary model, while the original A'^-body sim¬ 
ulation (not a real galaxy, though) continued to evolve substantially. 

The latter example also reminds about another conceptual 
problem: we always use a certain set of constraints in our methods, 
but there might exist other possible factors that we do not take into 
account, but that may distinguish between otherwise similar models 
having different internal structure and dynamics. An obvious exam- 
ple is the non-uniqueness of deprojection of a non-spherical galaxy 
(iGerhard & Binne^ll996h . or the insufficiency of using only the 
first two velocity moments to constrain the mass distribution, which 
could be all eviated by taking the f ull line-of-sight veloc ity distri¬ 
bution (e.g. iMerritt & Sahal1l993l : iGebhardt et alj|2()()3l . In Sec- 
tion l4.4l we only used the information about particle positions, not 
velocities, to construct a model which turned out to be quite close 
to the original one in terms of velocity profiles; nevertheless, had 
we imposed additional kinematic constraints, the model could have 
different properties. This also highlights another complication, spe¬ 
cific to the Schwarzschild method. Namely, it provides a solution 
that satisfies the given constraints in the best possible way, but 
such a solution may not be unique. In practi ce, this is circumvented 
by adding some regulariza tion procedures dRichstone & Tremaine! 
Il988t ICretton et alJll999ll . but a rigorously justified method for 
choosing one particular solution out of many is missing. 

Yet another issue is related to the difficulties in interpreting the 
outcome of Schwarzschild modelling in terms of convergence, il¬ 
lustrated by the model of Section l43l Our practical experience sug¬ 
gests that for a well-behaved model, with a large margin of feasibil¬ 
ity and the number of orbits exceeding the number of constraints by 
a large factor (> 10, i.e. an “underconstrained” model), regulariza¬ 
tion tends to create a narrowly peaked distribution of orbit weights. 
Conversely, if there are only a small number of orbits that get a high 
weight, this usually indicates that the required constraints are hard 
or impossible to satisfy exactly, implying that the chosen density 
profile might not be possible to realize with a self-consistent distri¬ 
bution function. Note, however, that a poor choice of initial condi¬ 
tions can also result in bad convergence, and that the inclusion of 
additional (e.g. kinematic) constraints may interfere with the ability 
of the solution to exactly match the di scretized density profile; for 
this reason, the Schwarzschild code of ivan den Bosch et alj j2008l) 
tolerates a few percent relative error in cell masses. 

There are many different methods for constructing equilib¬ 
rium models of disc galaxies, which we have reviewed in Sec¬ 
tion and tested some of them in Section ItTI but very few are 
suitable for non-axisymmetric (e.g. barred) galaxies, and no other 
implementation of Schwarzschild method is capable to deal with 
full three-dimensional models of multicomponent triaxial galaxies 
with figure rotation. Although our code in the current version is 
not intended for creating models based on observational data, it 
nevertheless presents a first step towards this task. Most models 
of external galaxies and even the Milky Way are created using the 
approximation of axisymmetry, which may introduce substantial 
biase s in recovered galaxy parameters, such as the mass-to-light 
ratio jThomas et al.ll200^ ; iLablanche et an i2012|) and the mass of 
the c e ntral supermassive black hole l lvan den Bosch & de Zeeuwl 
l201(]|;l0nkenetalj|2014l) . Moreover, the assumption of constant 
coefficient jSm of velocity anisotropy in the m eridional plane, co m- 
monly used in axisymmetric Jeans analysis dCappellarilbOOSh . is 
clearly violated for most models that we have considered. Thus 
having a more flexible modelling method will be increasingly im¬ 
portant for the current and forthcoming large galactic sur veys, such 
as ATLAS®° jCappellari. Emsellem & Krainovicll201 ih or CAL- 


IFA iSanchez et al.l20l3) . that yield a wealth of kinematic data for 
a large sample of galaxies. From the theoretical side, the interac¬ 
tion of non-axisymmetric disc s with triaxial haloes presents a num¬ 
ber of interesting effects (e.g.| B erentzen, Shlosman & JogedfioQ^ ; 
iMachado & Athanassoula! I^iol) . thus a flexible Schwarzschild 
method for creating such composite models is a valuable asset. 

5.2 Applications and implications of our Schwarzschild 
models 

5.2.1 Axisymmetric disc-bulge-halo models 

Simulat ions have s hown that galactic discs are prone to bar insta¬ 
bilities ( lHohllll97lh . with a growth rate which is function of many 
galaxy properties, such as halo, bulge, or gas m ass, velocity disper¬ 
sion in the disc etc. (see lAthanassoulall2013l . for a review). Thus 
such discs are not stationary. Yet in Section ItTI we showed that 
the SMILE software is capable of constructing axisymmetric multi- 
component models of disc galaxies. How are these two statements 
compatible? 

It is interesting to note that, although we did not set any kine¬ 
matic constraints, all the disc-bulge-halo models webuilt have high 
velocity dispersions in the disc, as can be seen in Fig. For all 
models, the minimum Q value is between 1.4 and 1.7, and occurs 
at radii between one and three disc scalelengths. Moreover, at radii 
less than one and more than three disc scalelengths, the Q values 
are mu ch larger (see Fig. [3]>. This ensur e s that the bar forms very 
slowly d Athanassoula & SellwoodllT986l ; lAthanassoulall2003h and 
the model is sufficiently close to stationarity for a Schwarzschild 
model to be possible. 

The fact that, in the absence of kinematical constraints, the 
Schwarzschild method finds a relatively hot disc as a solution, 
underlines the necessity of such constraints when modelling real 
galaxies and also suggests that it may not be possible to model cold 
discs using this method, unless the bar instability is slowed down 
sufficiently by the halo, bulge or gas. This leaves a sufficiently large 
range of axisymmetric disc galaxies to which the Schwarzschild 
method can be applied, including the lenticulars, the early-type spi¬ 
rals, and in general any hot discs. 

5.2.2 Rotating triaxial ellipsoids 

In Subsection 14.21 we built models for cuspy triaxial ellipsoids, 
which can be used to represent elliptical galaxies. While we have 
not attempted an extensive exploration of the relevant parameter 
space, our results indicate that, from the dynamical point of view, 
cuspy triaxial elliptical galaxies may exist in a broad range of pat¬ 
tern speeds. 

Observational studies indicate that elliptical galaxies fall into 
two classes - slow and fast rotato rs, with the former more likely to 
be triaxial temsellem et al.|[2007ll . In terms of the kinematic rota¬ 
tion parameter Ar introduced in that paper, our models cover the 
whole range from slow to fast rotators (Fig. lU. If we consider only 
models with no other constraint than the density distribution and 
the degree of figure rotation, then we find a rough general trend 
in the same sense as observations, i.e. more triaxial galaxies have 
smaller Ar values. If, however, we introduce further constraints, 
such as minimizing, or maximizing the global angular momentum, 
then the difference in the Ar value is much bigger than that due to 
a change of the pattern speed. 

Thus it becomes clear that further work on such models is 
necessary, e.g. to extend the axisymmetric Schwarzschild models 
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of lCappellari. Emsellem & BaconI (l2007h and answer whether it is 
possible to have strongly triaxial configurations with considerable 
figure rotation and large rotation parameters. 

5.2.3 Barred disc galaxy models 

As a next step, we attempted to use the Schwarzschild method 
to construct models of galaxies composed of a disc, a halo and a 
strong bar component. The potential and density of each of these 
components were given by analytic functions, as in orbital struc¬ 
ture studies. Although we made a large number of trials, we failed 
to construct a stable, self-consistent such model. More than one ex¬ 
planations can come to mind for this. 

One possibility is that the simple superposition of a number 
of rigid components may be far from a reasonable dynamical sys¬ 
tem. Indeed, all these components interact gravitationally and that 
could substantially modify their density distributions. For example, 
it is known that discs do not stay axisymmetric in the presence of 
bars, and form spirals, rings, density minima/maxima around some 
Lagrangian points, etc., while haloes will not stay spherically sym¬ 
metric at small radii if the potential of the thin discs dominates in 
the inner parts of the galaxy. 

A second alterna tive, which to our eyes is th e m ost likely one, 
is that , as discussed in lAthanassoula et al.l j20l4) an d iAthanassoulal 
( l2015h . the vertical mass distribution of Ferrers ellipsoids is unreal¬ 
istic. The real shape of bars is much more complex than a spheroid, 
having two components, one long and t hin (both horizontall y and 
vertically) and the other short and thick l lAthanassoulall200^ . The 
latter is also known as the boxy or peanut bulge, or the barlens. 
Since the orbits in the bar potentia l are able reproduce this complex 
bar s hape (e.g. [Pfennigea Il984al : IPatsis. Skokos & Athanassoulal 
l2002h . they may not be able to reproduce the simple ellipsoidal 
shape of the Ferrers bars as well. A further argument in favour of 
this alternative is the fact that it has been possible to make two- 
dimensional mode ls of barred galaxy potentials, as shown both by 
IPfennige^ (Il984tit) and in Section l4^ Moreover, in Section l4!4l we 
have been able to make a Schwarzschild model for barred galaxies 
having a realistic vertical distribution. 

5.2.4 Barred N-body galaxy models 

We also attempted to build Schwarzschild models for two W-body 
snapshots taken at the beginning of the secular evolution phase 
of the bar and using only the densities as a constraint. One of 
the two cases was specifically chosen from a simulation which 
showed hardly any bar evolution, while the other, on the contrary, 
was chosen from a simulation with a very strong evolution. For 
the former we were able to build a self-consistent model, which, 
when evolved, followed the evolution of the snapshot. For the lat¬ 
ter, however, although we could build a self-consistent model, its 
evolution was considerably less than that of the simulation. This 
presumably means that, of all the velocity distributions compatible 
with the density distribution in the simulation, in absence of kine- 
matical constraints, the Schwarzschild model found the one which 
was nearest to stationarity (see also the corresponding discussion 
in Section [5.2.1b . 

This may set some limits to the applicability of the 
Schwarzschild method to barred galaxies, if barred galaxies are 
found to have non - neglig ible evolution. Indeed, as shown by 
lAthanassoulal 1 I 2 OO 2 L 1200311 . the resonances in the halo absorb a 
large fraction of the angular momentum emitted from the bar re¬ 
gion. Since no kinematic constraints can be set on the haloes, the 


Schwarzschild method may tend to provide self-consistent models 
fulfilling all the density constraints as well as some kinematic con¬ 
straints for the disc, but whose evolution may be totally different 
from that of the real galaxy. Thus more work is necessary here be¬ 
fore the Schwarzschild method can be massively applied to barred 
galaxies, including our own Milky Way. 


6 CONCLUSIONS 


We have presented an implementation of the Schwarzschild or¬ 
bit superposition method for creating equilibrium dynamical mod¬ 
els of non-spherical stellar systems. Thi s is a further de velopment 
of the publicly available SMILE code dVasilievI l2013h . with the 
major new and improved features primarily targeting disc galax¬ 
ies: a new general-purpose potential approximation in terms of 
two-dimensional spline-interpolated Fourier expansion in cylin¬ 
drical coordinates, which is accurate even for very flattened sys¬ 
tems, full support for multicomponent systems and for figure ro¬ 
tation. The software is suitable not only for constructing self- 
consistent models with given properties, but also for analysing the 
orbital structure of a given potential, ei ther specified analytically or 
taken from an A^-body simu lation (e.g. lHarsoula & Kalapotharakoa 
I 2 OO 9 I : lAthanassoulal I 2 OI 2 I) . This opens up new possibilities for 
studying the properties of orbits relevant for the formation of bars, 
rings, spiral arms, etc., which traditionally have been explored by 
integrating the motion of test particles in a given analytic potential 

e tos. Patsis & Athanassoula l 2002l : [^mero-G6mez et alj201ll : 

rri. Antoia & Flelmil I 2 OI 3 I) . and could now be improved by 
taking these particles from a self-consistent model while retain¬ 
ing the smooth potential for the orbit calculation. One may even 
use several snapshots from an A-body simulation to construct a 
time-dependent smooth potential approxi mation and use it for or¬ 
bit integration l lManos & Machadal2014ll . To facilitate these ap¬ 
plications, the module for computing potential and forces is pre¬ 
sented as a separate library SMILEPOT, with C and PYTHON inter- 
faces and bind ings to oth er stellar-dynamical software tools: NEMO 
(Teuben 19951) . AMUSE jPortegies Zwart et al.ll2013l) and GALPY 
lBovvll2015l). 


We have confirmed the possibility of constructing steady- 
state models of disc galaxies, both axisymmetric and barred, with 
the Schwarzschild method. We compared several existing codes 
for creating composite axisymmetric disc-bulge-halo models and 
found that SMILE performs very well in this task, producing mod¬ 
els that are closest to equilibrium. While elliptical galaxies are not 
the main subject of this paper, we also considered models resem¬ 
bling triaxial elliptical galaxies with density cusps and figure rota¬ 
tion. We have shown that they may exist in a wide range of pattern 
speeds and angular momenta, and are stable over many rotation 
periods. On the other hand, a barred disc model with an analyti¬ 
cally defined potential, extensively used for orbit analysis in previ¬ 
ous studies, turned out to be impossible to realize self-consistently, 
and the closest approximation to it was quite unstable dynamically. 
Nevertheless, there exist models of barred disc galaxies embedded 
in nearly-spherical haloes, that appear to be in almost steady state. 
We have constructed two such models using density profiles ex¬ 
tracted from A-body simulations, utilizing only information about 
particle positions. These models were shown to be reasonably sta¬ 
tionary, and their kinematical properties were similar to the original 
A-body models if the pattern speed matched these original simu¬ 
lations. Interestingly, one of the original A-body models demon¬ 
strated considerable evolution (slowdown and growth of the bar). 
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while the corresponding Schwarzschild model evolved much less. 
We explore the reasons for this difference in a future work. 
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file). Let the density be represented by a Fourier series 

'^max 

p{R,z,(j>) = ^ pm{R,z) exp{im(j)), 

m=0 

so that 

r2-K 

pm{R, z) = {2n)~^ / p{R,z,(l))exp{-im(j)). 

Jo 

The potential generated by this density is represented as 

^max 

${R,z,cf))= ^ $m{R,z) exp{im(l)), 

m=0 

with each term in the expansion 

/ + 00 POO 

dz dR'2nR' pm{R\z')E, 

-oo J 0 

PCG 

E{R, z, R', z') = / dk Jm{kR)Jm{kR') exp{—k\z — z'\). 

Jo 

The last integral can be evaluate d malytically 

dOradshtevn & Rvzhikl 1 196.1 ICohl & Tohli'^ Il999h . using 
the Legendre function of the second kind Q\ 


r—-—Qm—1/2 ( 

tVrW V 


R^ + R'^ -hiz- z'f 
2RR’ 


if R, R' > 0, 


H ^{R^ -F + (z - z'f)-^^^ Smo if i? = 0 or i?' = 0. 


For a discrete point mass set, the potential coefficients are 
computed as 

N 

^m{R,z) = —G'^mpE(R, z, Rp, Zp) expimclip. 

p=i 

The spline potential expansion in cylindrical coordinates 
is initialized by computing the coefficients ^rn{R, z),m = 
0 ... mmax either from a given analytic density profile, or from 
an A'^-body snapshot, on a two-dimensional grid in R, z cov¬ 
ering a rectangle R < i?max, 1^1 < Zmax- To improve reso¬ 
lution at small radii, we use logarithmic scaling of coordinates 
(R = ln(l -F R/Ro), and similar for z), and to reduce the ap¬ 
proximation error at large radii, we scale the potential by a factor 
fR^o + R^ + z^ so that for large R, z the scaled coefficient for 
m = 0 tends to a constant. Flere Rp = —GMtotai/3>o(0, 0) is a 
characteristic radial scale of the system. The evaluation of poten¬ 
tial, forces and density at an arbitrary point within the box is done 
by computing interpolating splines or their derivatives. Outside the 
box, density is assumed to be zero, and the potential and forces are 
computed using a quadrupole approximation to the mass distribu¬ 
tion. 

Fig. lAll demonstrates that for strongly flattened models, the 
cylindrical spline expansion is superior to the spherical-harmonic 
expansion, although if the potential is initialized from an A^-body 
snapshot, the accuracy is mainly limited by discreteness noise. 


APPENDIX A: POTENTIAL EXPANSION IN 
CYLINDRICAL COORDINATES 


To evaluate the potential at an arbitrary point in cylindrical 
coordinates R, z, rj) , we may use the following algorithm (see 
Binnev & Tremaina 120081 section 2.6.2 for razor-thin discs and 


Cuddefordll993 for axisymmetric discs with arbitrary vertical pro- 
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Figure Al. Mean square error in potential (left), force (centre), and density (right) approximation, as functions of radius. 

Top row: axisymmetric Miyamoto-Nagai profile with A = 3, S = 1, using spherical-harmonic expansion (dashed lines) with the number of angular terms 
ranging from 10 to 50, and using spline interpolation on a grid in the meridional plane (dash-dotted lines), with different number of nodes (grid extends to 
R = 2000, 2 : = 100). Clearly, the latter provides much better accuracy than the spherical-hai'monic expansion for the strongly flattened density profile, which 
converges very slowly. 

Middle row: the same density profile sampled by = 10^ points. In this case, the accuracy of cylindrical spline is limited by the discreteness noise at 
all radii (increasing the grid size does not make it smaller), while for the spherical-harmonic expansion the number of angular terms still limits the accuracy 
at large radii, where the density profile is very flattened. Note that both smooth potential expansions still approximate the target model better than a direct 
evaluation of potential and density from the A^-body snapshot. 

Bottom row: strongly triaxial Ferrer s model with axes x : y : z = ^ \ l .b : 0.6. Here, we compared the spherical-harmonic expansion with the radial basis set 
from spherical Bessel functions (as in i Allen. Palmer & Papaloizoill99(ll) with 40 radial terms, and the cylindrical spline with the grid size = 20, Nz = 15; 
in both cases, the number of angular terms was 10. Again the latter was found to approximate the target model better than the former. Left half of each panel 
compares the approximations initialized from a smooth density profile, while the right half does the same for a. N = 10^ particle snapshot; again in the latter 
case the discreteness noise limits the accuracy. 
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